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ABSTRACT 

We use the weighted integral form of spherical Bessel functions, and introduce a new 
analytical set of complete and biorthogonal potential-density basis functions. The 
potential and density functions of the new set have finite central values and they fall 
off, respectively, similar to r^(^+'^ and r^^'^+^'> at large radii where I is the latitudinal 
quantum number of spherical harmonics. The lowest order term associated with I = 
is the perfect sphere of de Zeeuw. Our basis functions are intrinsically suitable for the 
modeling of three dimensional, soft-centred stellar systems and they complement the 
basis sets of Glutton-Brock, Hernquist & Ostriker and Zhao. We test the performance 
of our functions by expanding the density and potential profiles of some spherical and 
oblate galaxy models. 

Key words: celestial mechanics, stellar dynamics - galaxies: kinematics and dynam- 
ics - methods: analytical - methods: numerical 



1 INTRODUCTION 

Solving Poisson's equation is an im portant step in the 
, study of self-gravitating stellar systems (|Binnev fc Tremaind 
120081 ). Expanding the density distribution and its con- 
jugate potential field in terms of a complete basis 
set is one of the most efficient methods that inves- 



onal potential-density (PD) basis sets have been developed 
by Clutton-Brock (1973, hereafter CB73) , Hern quist & Os- 
triker (1992, hereafter H092) and IZhad (|l996l '). CB73 and 
H092 set the lowest order terms of their b asis functions 
to th e lPlurnmej ImJ) and iHernguistl l|l990l ) models while 
IZhad I 19961 ) uses an Q-model with the density 



tions (iFridman & Polvachenko 


Il984l: iHernauist & Ostrikeij 


19921: Earn & SellwoodI 1995 


; Meza & Zamorand 19971: 


Weinberg & Katd l2007l: iBuvle et all |2007|) and the first- 



p(r) 



C 



(1) 



order stability analys is of both fla t (iKalnaisI 119771 : 
IPichon fc Cannonlll997l : Ijalali fc Hunter 2005: Jalah 2007) 
and three dimensional galaxies l(Saha.l991. : .Weinberg.! 1991) . 
Consequently, the success of those studies highly depends 
on the choice of basis set. Desirable potential and density 
basis functions should be biorthogonal and converge rapidly 
in order to decrease the computational noise and cost. Nev- 
ertheless, finding a suitable basis set is not an easy task 
and only few analytical basis sets have been found for three 
dimensional stellar systems. 

For stellar systems of finite size, spherical Bessel 
functions are the classical biorthogonal eigenfunctions 
of the Laplace operator and they have been used 
in the stability analysis of certain spherica l galaxies 
IIFridman fc Polvachenkolll984l : lAllen et al.lll990l : IWeinbergj 
Il99ll ). For galaxy models of infinite extent, three biorthog- 
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where C is a constant parameter. The lowest order term 
of a PD s et does not necessarily need to be spherical 
dSverl Il995h . but an orthonormalization using the stan- 
dard Gram-Schmidt procedure must be adopted (|Sahalll9^ : 
iRobiin fc Earn|[l996l ) to guarantee the completeness of the 
set. 

Apart from the quoted analytic basis functio ns, numer- 
ically generated sets have also become available. IWeinberS 
(|l999t ) assumed the form of the lowest order basis func- 
tions and numerically solved the Strum-Liouville equation 
to obtain biorthogonal basis functions of higher orders. De- 
spite this worthwhile contribution, the propagation of com- 
putational noise during the application of numerical basis 
functio ns has become problematic in recent A''-body experi- 
ments (jKalapotharakos et al.l[2008l ). which justify the ongo- 
ing search for new analytical basis functions. 

In this paper we introduce a new analytical set of 
biorthogonal PD basis functions whose potential and den- 
sity components have finite central values, fall off similar to 
H092 functions a s r — > oo, and the ir lowest order term is the 
perfect sphere of Ide Zeeuwl l|l985h . We derive and evaluate 
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the weighted integral forms of spherical Bessel functions in 
311 and obtain the radial basis functions in terms of associ- 
ated Legendre functions. In !j3l we use the new basis set and 
generate the series representations of certain spherical and 
oblate galaxy models. We end the paper with concluding 
remarks. 



2 POTENTIAL-DENSITY PAIRS 

We define r = (r, 6, (p) as the position vector expressed in 
terms of usual spherical coordinates, with r, 6 and </> be- 
ing the radial distance from the origin, co-latitude and az- 
imuthal angle, respectively. We also assume that the mean- 
field potential and density functions of a stellar system ad- 
mit the following expansions 



nlm 

p(r) = AiimPnim(r) 



(2a) 
(2b) 



The basis functions $nim(r) and p,iim(r) satisfy Poisson's 
equation 

V^$„i™(r) = 47rG'p„,„(r), (3) 

where n and I are the radial and latitudinal quantum num- 
bers corresponding to r and 6, respectively, and m is the 
azimuthal Fourier number associated with (j>. G is the uni- 
versal constant of gravitation. We proceed with a case that 
Pnim is proportional to ^„im- This reduces Poisson's equa- 
tion to the eigenvalue problem 



(4) 



that involves the Laplace operator V'^ and a constant pa- 
rameter k. Since the Laplace operator is Hermitian, its as- 
sociated eigenfunctions form a complete biorthogonal basis 
set. The coefficients Pnim and D„im thus become identical. 
The spherical harmonics Yim{9,4') and the spherical Bessel 
functions jnikr) are the classical solutions of (|4]). 

While Yim.{6, <t>) show an acceptable performance in the 
expansion of physical quantities in terms of angle variables, 
Bessel functions do not look lik e galactic profile s and can not 
generate efficie nt expansions llWeinbereflT999l ). We extend 
the method of IClutton-Brockl (|l97^ to three dimensional 
systems and express the eigenfunctions as 



where 



$,i!m(r) = -Ylra{0,(t))^nl{r), 
Pnlra{r) = Ylra{0 , (}>) p„l{r) , 



^ni{r) = / ji{kr)g„i{k)dk, 



(5a) 
(5b) 



(6a) 



Pni{r) 



AttG 



ji{kr)g„i{k)k dk. (6b) 



The functions gni(k) (n, / = 0, 1, 2, • • • ) are to-be-determined 
functions that we require to satisfy the biorthogonality con- 
dition 

/ $,i!m(r)[p„'i'm'(r)]*dr = (7) 
Here, the asterisk denotes complex conjugation and 5^/ is 



the Kronecker delta. Substituting from ((5)1 and ((Sjl in ((Tjl 
and using the identity 

fl r2TT 

d(cose) / d4,Y^{e,4,)Yi%^.{e,^) = 5w5rr,m', (8) 



-1 Jo 
the orthogonality condition ((7|) reduces to 

roc roo 

AnG J, 9ni{k)dk J g^n {k') k'^dk' 

/oo 
ji{kr)ji [k'r) r^dr = InimS„„' 



(9) 



The innermost integral on the left-hand side of S is eval- 
uated according to the Fourier-Bessel theorem ifuginciusl 
11972) as 

ji{kr)ji{k'rydT=^5{k' ~k), (10) 

with 5{k' — k) being the Dirac delta function. Substituting 
(fTOl) in (O leads to 



8G 



gni{k)g„'i{k)dk = Ini5„„i = I„im5„„'- (11) 



This condition requires gni{k) to be any orthogonal set of 
functions over the semi-infinite fc-domain. Our special choice 
is gni{k) = k'L^{2k)e~'' where L^{k) axe the associated 
Laguerre polynomials that obey the following orthogonality 
relation 

e-*'k^Ll(k)Ll,{k)dk^^^^^5,,,. (12) 

Consequently, the constant parameters on the right-hand 
side of equation (lllf) become 

- _ (n + 20! 
^ G22i+4n! ' 
and our radial basis functions read 

V'ni(r) = / jiikr)L^J{2k)e-''k'dk, 
Jo 



(13) 



(14a) 



Pni(r) = — ji{kr)Lt:i2k)e-'k'+'dk. (14b) 

The integrals in (I14p converge rapidly, which makes 
their evaluation a straightforward task by numerical meth- 
ods. Ifowever, closed-form analytical expressions can also be 
derived for tpni{r) and p„i{r) as we explain below. We utilise 
the series form of the Laguerre functions 



L^'(2fc)=^(-1) 



n + 2l\ {2k) 



(15) 



and express ji{kr) in terms of Bessel functions to rewrite 
(ll4bD in the form 



1 " 



, /n + 21 



X / J,+ i{kr)e "k 



TT_1 

2r i\ 

k,l+i+3/2 



dk. 



(16) 



Carrying out a change of independent variable as kr u, 
transforms equation p6p to 



pnl( 



i f n + 2l\ jTxr 



X / J,, i(M)e-"/''u'+'+^/''du. (17) 
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Figure 1. Left panels display the radial parts of tlie potential basis functions, 4'nl{^)} for ^ = 0,1,2,3,4, and right panels show their 
conjugate density functions pni{r)- Top and bottom panels correspond to i = and 1 = 2, respectively. All functions have been normalised 
to their maximum values. 



T he intepjral in (|17ll can b e calcu lated using equation (6.621) 
in iGradshtevn fc Ryzhi^ ||2000[). Defining ^ = l/VT+r^, 
vi = I + i + 3, U2 = Z + i + 3/2 and ^ = —1/2 — I, we obtain 



which start from 



■E 



(2Z + i + 2)(2« + i + l)(n + 20! 



4V27rG^ i\{n-i)\ 
^2\-i/4 



{-2yC'{^-e) ' P!^AO, (18) 



where Pj^ [x) are associated Legendre functions. Following a 
similar procedure, one can show that 

' i=0 

(19) 

where vg, = I -\- i -\- 1 and V4, = I + i — 1/2. The associated 
Legendre functions can be determined through the recursive 
relations (jCradshtevn fc Rvzhil3l2OO0h 

- M + i)J'^Vi(e) + ('^ + = (2^^ + mP.'iO, 

(20a) 

J'^i(e) - ^'."+1(6 = (2;.+ 1)(1 -C')i/^Pr'(0, (20b) 



-1/2 



(0 = 



2 sin [(;/ + 1/2) /?] 
TTsin/? (i^ + l/2) 



cos/3 = C. (21) 



The lowest order members of our PD family are 
1 



^oo('') ~ arctanr, 

r 

1 1 

Poo(r) 



(22a) 

2^( TT^' 

which define the perfect sphere of lde Zeeuwl (Il985l ). We have 
therefore found a biorthogonal basis set that is distinct from 
CB73, H092 and Zhao's (1996) functions. Moreover, from 
(|21|) and the recursive relations (|2Up . we deduce that the 
functions P^^(C) and P^^i^) behave, respectively, similar to 
r° and in the limit of r — > 00. It can thus be verified 
that ipni{r) ~ r"(^+'' and pniir) ~ r"(*+') hold at large 
radii. The potential functions of CB73, H092 and ours have 
finite central values and they fall off similar to r^*^^^'' at 
large radii. Our density functions are analytic at the galac- 
tic centre as are the functions of CB73, but they behave 
like H092 functions in the limit of r — > cx). The best per- 
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formance of our basis set is thus expected in soft-centred 
systems whose outer density profiles are similar to r"*. 

In Figure [1] we have displayed several members of our 
basis functions for I = Q,2. At the centre, both the potential 
and density functions have finite, non-zero values for / = 0, 
and they vanish there for I ^ 0. The expected yet interesting 
property of V'ni (r) and pni (r) is their oscillatory nature. The 
number of peaks of our functions (in the radial direction) is 
equal to n -I- 1. Our numerical experiments show that the 
series built by oscillatory functions have a faster and more 
accurate mean-convergence compared to functions that do 
not share this feature. 

Our functions have a length scale that has been set to 
unity so far. In general, changing the length scale is neces- 
sary to reconstruct galaxies of different core radii. A scal- 
ing parameter ro can be easily introduced t o our formula- 
tion through replacing gni{k) with gni{kro) (|Clutton-Brockl 
Il972[ ). This implies the following transformations 

ro"'$n!m(r/ro,6',0), (23a) 
(r/ro,e,<^), (23b) 
(23c) 



10" 



Pnlm {r, 0, I 



3 RECONSTRUCTION OF MODEL GALAXIES 

Bi-orthogonal basis functions, similar to ours, have the ad- 
vantage that the coefficients Pnim = Dnim in ((2]) can be 
determined using either the potential $(r) or the density 
p(r) through the following formulae 



1 

Inl 



p(r)['l>„i,n(r)]*dr 



=- I 'I'(r)[p„i„(r)]*dr, 

J-nl 



(24) 



where we have used the orthogonality conditions (O and 
In what follows, we examine the performance of our 
basis functions by the series reconstruction of the density 
profiles and potential fields of some model galaxies. 



3.1 Spherical Models 

We followed the standard procedure of using spherical har- 
monics for the expansions of physical quantities in terms of 
angular variables, and introduced a new set of radial basis 
functions. So we need to examine the performance of our ra- 
dial set by reproducing some spherical models. As case stud- 
ies, we choose the isochron e and Plummer models of total 
mass M and length scale b (|Binnev fc Tremaindliooil ). Our 
basis functions have finite values at the centre and it would 
be interesting to learn whether they are suitable for the re- 
construction of models with central density cusps. For doing 
so, we also analyse the performance of our basis fun ctions 
by applying them to Dehnen's 7- models l|DehnerJ 19931 ). The 
density profiles of Dehnen's models diverge similar to r~~' in 
central regions and fall off proportional to r"'* at large radii. 
Dehnen's models also have a length scale h. The model with 
7 = has an intrinsic core at the centre and for 7 = 3/2 a 
centra l cus p with an in termediate slope between iHernguistl 
l| 19901 ) and lJaffd l|l983D models is created. In our study, we 
choose two models with 7 = 1/2 and 7 = 3/2. 



1 

C5 




Figure 2. The absolute magnitudes |DnOo| of the expansion coef- 
ficients versus the radial quantum number n for several spherical 
models. 



We have set ra = 1 and used equation (|24p to com- 
pute the coefficients of expansion D„im for the Plum- 
mer, isochrone and Dehnen models. The parameters of the 
isochrone model have been set to GM = 1 and b = 0.5. 
For other models we have used GM = 1 and b — 1. Our re- 
sults are displayed in Figure (2] which shows how \Dnoo\ vary 
versus n. We note that all coefficients with Z, m 7^ vanish 
because of spherical symmetry. It is evident that |-D„oo| de- 
crease several orders of magnitude by including more terms 
in the series expansions. Although for 7=1/2 the coeffi- 
cients of Dehnen's model fall off similar to other soft-centred 
models, they decay mildly for 7 = 3/2. This shows very 
slow and unfavourable convergence of our series expansion 
in steeper cusps as is expected. 

Having expansion coefficients, the original model can 
be constructed using Denoting the original PD pair by 
[$o('"), po{r)] and their series representations by [^{r), p{r)], 
we compute the relative errors Ep—{p—po)/po and £# = ($ — 
$o)/4'o, and their absolute magnitudes 5p = \Ep\ and (5# = 
\E^\ to measure the performance of the basis set. We have 
used the first 10 radial basis elements (rimax = 9) to compute 
[$(r), p(r)]. The results are shown in Figure[3] It is seen that 
5$ is below 2% in all parts of the Plummer, isochrone and 
Dehnen's 7 = 1/2 models, and also for r > 0.1 in Dehnen's 
7 = 3/2 model. The reason is the similarity of tp„o{r) defined 
in (|19p to the potential profiles of the chosen models. The 
large error magnitude near the centre of Dehnen's 7 = 3/2 
model is due to its sharper density cusp that prohibits a 
simultaneous convergence of the density and potential series. 

The reconstruction of po{r), however, has not been suc- 
cessful in Dehnen's 7 = 3/2 model because of its sharper 
cusp. Large values of 5p are also observed in the central part 
of Dehnen's 7 = 1/2 model (due to its cuspy nature), and 
at large radii of the Plummer model due to its rapid den- 
sity fall-off, which is steeper than our p„o(r) ~ r""*. The 
isochrone model is the only case that has been reproduced 
with a reliable accuracy in all parts of the galaxy. In fact, 
the isochrone model shares two basic features of our new 
basis set: (i) It has a soft core, (ii) Its outer potential and 
density profiles decay, respectively, similar to and r~* 
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Figure 3. The relative errors 5$ (tiiin solid lines) and 5p (thin dash-dotted lines) in the series reconstruction of the potential |<l?o| (thick 
solid lines) and its associated density po (thick dash-dotted lines). The basis functions are those of equations I II8I I and II19I I. Panels a, b, 
c and d correspond, respectively, to the isochrone, Plummer, Dehnen's 7 = 1/2 and Dehnen's 7 = 3/2 models. The functions <I>(r) and 
p(r) have been computed by taking the first 10 basis functions (rimax = 9)- In each error curve, there are rimax + 1 sharp minima that 
correspond to the locations of exact match between the original functions and their series representations. The scattered squares show 
the magnitude of 5p calculated from equation 1291 1. The central values of $0 and po have been normalised to some arbitrary numbers 
just for the purpose of visualising their profiles against the error curves. 



as do the envelopes of the functions ?/'„o('") and p„o{r). In 
Figure [Sja for rimax = 9, the magnitude of Sp is less than 
1% over the range < r < 20 and it saturates at a level 
of Sp ~ 5% for r > 20. By increasing rimax to 14, both 5p 
and S<s> remain smaller than 1% over the range < r < 100. 
This result can also be deduced from Figure [2] that shows a 
monotonic decay for \DnOo\ versus n. 

In general, the density error Sp is larger than 5$. We 
explain this by calculating Ep in terms of Ei> and its 
derivatives. The original potential and density functions sat- 
isfy Poisson's equation, and since our basis functions are 
biorthogonal, the relation V^cl? — AnGp also holds between 
the expanded quantities. We can therefore write 



V' ($-$0) =47rG (p-po): 



(25) 



which is divided by po to obtain 

— ($ - $0) = -Lv^ (^oSs) = ^-wGEp. (26) 
Po Po 

For the Laplace operator with spherical symmetry, equation 



pS)) leads to 
1 



po 



$oV^£4> + E^V'^o + 2—^—^ 
ar ar 



= 4ttGEp. (27) 



Substituting 47rGpo for V^^o in (p7)) . yields 



^V^E^ + 4nGE, + ^^^ 
Po po ar 



dr 



4nGEp 



(28) 



We are interested in the local extrema of E<i, . There are n^ax 
number of such points whose existence is deduced from the 
oscillatory nature of basis functions. The derivative dE^/dr 
vanishes at the extrema of i5# and equation (|28p reads 



$0 



1 d^E<i 



E^ AnGpo E^ dr^ ' ^ ' 

This is a useful relation that gives a credible estimate of 
Ep/Eis> based on the quotient $o/po and the curvature of 
For each model, we have independently computed Sp 
from (|29|) and have plotted the results (scattered squares in 
Figure [3)l against the numerical graph of Sp obtained from 
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Figure 4. Variations of 5p {top panels) and 5^ (bottom panels) for the isochrone {left panels) and Dehnen's 7 = 1/2 {right panels) 
models. Solid, dotted and dash-dotted lines respectively correspond to expansions by our, CB73 and H092 basis functions. The length 
scale of the isochrone model is 6 = 0.5 and that of Dehnen's 7 = 1/2 model is 6 = 1. For both models we have set GM = 1 and all basis 
functions have the length scale of ro = 1 . 



the series expansion of po{r). There is a close agreement 
between the results of two methods, confirming the fact that 
the drift 5p — 5$ is independent of the choice of basis set and 
it persists in any series solution of Poisson's equation. 

Neither the isochrone nor Dehnen's 7 = 1/2 models 
match the zeroth order terms of CB73, H092 and our basis 
functions. Therefore, the performance of these basis sets can 
be fairly compared by expanding the isochrone and Dehnen's 
7 = 1/2 models (Figure |4]). It is seen that CB73 functions 
have a poor performance in reproducing both models. Our 
functions have performed better than H092 functions for 
r < 1 in the isochrone model. Nevertheless, H092 functions 
have resulted in the lowest magnitudes of 5$ and 5p for r > 
1 in the isochrone model, and for r > 0.005 in Dehnen's 
7 = 1/2 model. Our results show that the envelopes of basis 
functions must follow the radial profiles of both the density 
and potential functions of a spherical stellar system to assure 
a reliable expansion. Dehnen's shallow density cusp cannot 
be reproduced even by cuspy set of H092 (see Figure |4f)) 
because the central envelope of H092's density functions is 
proportional to while Dehnen's density profile diverges 

r.-l/2 



3.2 Oblate Galaxy Models 

The modeling of oblate galaxy models is a bigger challenge 
because the series of radial basis functions must converge 
together with spherical harmonics. It is therefore hard 
to predict how the combination of radial and angular 
functions will behave. As our case s tudies of spheroidal 
galaxy models, we choose an oblate iKuzmin fc Kutuzovl 
l|l962f) model and a perfect spheroid (|de Zeeuwl 119851 )7 
and reproduce their density distributions using the series 
of CB73, H092 and our new biorthogonal sets. In the 
spherical limit, the Kuzmin-Kutuzov model reduces to 
Henon's (1959) isochrone, and the perfect spheroid becomes 
the perfect sphere, which is the lowest order term of our 



new basis set. Since our functions showed slow convergence 
for Dehnen's spherical models near the centre (see Figures 
[2] an d[3l), we did not extend our analysis to their flattened 
(Dc hnen fc GerhardI 1 1994 ) counterparts. Moreover, we 
proved in i^S.ll that the potential expansions are always 
more accurate than the density ones. This applies to oblate 
models as well, and therefore, we confine ourselves to 
computing Sp. 

Defining = a^c^ + c^B? + a^z^, the density functions 
of the Kuzmin-Kutu zov and perfect spheroidal models are 
respe ctively given by l|Deionghe fc de Zeeuwl 198^ : Ide Zeeuwl 
119851 ) 



pKKiR, 



Mc^ (a^ - c^)R'^ +a* + 2u^ + Sa^M 

47r u3 (i?2 j^z'^ + a^ + c^ + luf^ ' 

-2 



PPS(-R,2:) = 



M 



1 + + 



1-6= 



(30a) 



(30b) 



7^2^/^^ 

where R = rsinS and z = rcos9, and z is the symme- 
try axis. The parameter e is the fiattening of the perfect 
spheroidal model. The Kuzmin-Kutuzov model has equipo- 
tential surfaces of the axis ratio \fcfa near the centre, and 
we choose its length scale so that a 4- c = 1. Here again, M 
is the total mass of the galaxy. 

The isocontours of the original and expanded density 
functions are displayed in Figure [S] for a perfect spheroidal 
model of e = 0.5 and for a Kuzmin-Kutuzov model of 
c/a = 0.5. We have set ro = 1 and nmax = 9, and used 
Zmax = 8 for the Kuzmin-Kutuzov model and Zmax = 4 for 
the perfect spheroid, respectively. The maximum deviation 
from the original model occurs near the J?-ax;is because of 
the slow convergence of spherical harmonics a& Q ^ -n 12. 
Therefore, we have shown in Figure |6] the variation of 5p 
versus R in the equatorial plane. Note that the existence 
of a symmetry axis implies Dnim ~ for odd latitudinal 
quantum numbers and for m 7^ 0. 

Our experiments show that by increasing Zmax the den- 
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0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 

R R 



Figure 5. Left panel: Density isocontours of a perfect spheroid. Solid lines correspond to the exact model density and short-dashed, 
long-dashed and dotted curves are associated with density expansions using our, CB73 and H092 basis functions, respectively. The 
ellipticity of the model is e = 0.5. Right panel: Same as the left panel but for an oblate Kuzmin-Kutuzov model with c/a = 0.5. In both 
figures, the levels of isocontours indicate the fraction of maximum density. 




Figure 6. Radial variation of the relative error Sp (oscillatory curves) against the model density on the equatorial plane (z = 0) of the 
same prefect spheroid (panel a) and Kuzmin-Kutuzov (panel b) models of Figure [5] Dashed, dotted and thick solid curves correspond, 
respectively, to CB73, H092 and our basis functions. 



sity expansion near the equatorial plane is improved. It is 
evident that H092 functions have failed in reproducing the 
finite central densities of both models but they have best 
fitted the outer parts. For R > 0.2, the error indicator Sp is 
smaller for H092 functions than CB73 ones by almost one 
order of magnitude, and that of our new functions lies be- 
tween them. Nonetheless, only our functions result in very 
small error level of < 1% for R < 0.2 in both models. This 
shows that our new basis set is the most trusted tool for 
modeling all parts of cored, oblate galaxies whose outer po- 
tential and density profiles fall off similar to r^^ and r~^, 
respectively. We note that the magnitude of Sp rises substan- 
tially and then saturates beyond the radial distance R ~ 10 
where the density has fallen to 0.1% of its central value. This 
property is shared by all tested basis sets. It is by increasing 



the number of radial basis functions (rzmax) together with 
the precision of computations that the error magnitude is 
suppressed at large radii. 

It is helpful to compare our results with iRobiin &: EarnI 

1 19961 ) who have designed a set of basis functions for the per- 
fect spheroidal models. Their functions have been orthonor- 
malised using Gram-Schmidt procedure. For iimax ~ 9 and 
'max ~ 4 that match the number of series terms in our setup, 
they reported a maximum error of Sp ~ 30% in the density 
expansion for a perfect spheroid of ellipticity e — 0.5 and 
inside the domain < R, z < 2. In the same region, our 
density expansion leads to a maximum error of Sp « 7%, 
which is notably small. 
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4 CONCLUSION 

The lack of suitable PD basis sets is a serious problem in 
dynamical studies that solve Poisson's equation using series 
expansions. For three dimensional stellar systems only few 
analytic basis sets have been found and most researchers 
have tailored numerical functions to cope with their spe- 
cific problems. In this paper we generalised Glutton-Brock's 
(1972) idea to three dimensional systems and introduced a 
new set of basis functions, which have the useful property of 
biorthogonality. Our functions complement the CB73, H092 
and Zhao's (1996) basis sets because neither of them exhibits 
the following properties together: (i) A finite central density, 
(ii) An outer density fall-off similar to r^*. For instance, the 
integrable models of de Zeeuw (1985) and their perturbed 
states, can be efficiently expanded by our basis functions. 
Thus, we get one step closer to the stability analys is of el- 
liptica l gal axies whose potentials are of Stackel form. iRobiinI 
l|l995ll and lSellwood fc Valluril l|l997h investigated the insta- 
bilities of some spheroidal galaxy models but calculating the 
eigenspectra of more general triaxial systems remains as a 
big challenge. 

Our functions were derived in terms of elementary ratio- 
nal, and associated Legendre functions for which recursive 
formulae are available. We carried out a mathematical error 
analysis and then compared its results by numerical experi- 
ments to show that density expansions converge slower than 
potential ones. By expanding several spherical and oblate 
galaxy models, we showed that an improper choice of basis 
functions can contribute potentially dangerous errors to dy- 
namical studies. Not only the nature of the galactic centre 
(cuspy or cored) is an important factor for the selection of 
basis functions, the outer density and potential profiles also 
matter. Neither our new set, nor other basis functions cited 
in this paper, are suitable for the modeling of cuspy dark 
matter halos whose density profiles decay outward like . 
It is possible to find basis sets compatible with such systems, 
but that will require other choices of the weighting functions 
gniif) that must be orthogonal over the r-domain in three 
dimensions. It is remarked that we had set the length scale 
of our basis functions to ro = 1 in all of our case studies, 
but there is always an optimum value of ro that gives the 
best fit. For example, the isochrone model is best fitted by 
setting ro = 26. We therefore recommend an optimal search 
for finding the best minimiser of 5$ . 
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